Este reporte explora los resultados del monitoreo participativo del estado de salud de árboles de oyamel en el Parque Nacional Desierto de los Leones y sus zonas de influencia en Santa Rosa Xochiac.

Los datos corresponden al resultado del muestreo participativo realizado por 12 brigadistas de Santa Rosa Xochiac, utilizando kobo-conabio como parte del proyecto 308488 Monitoreo y manejo para la conservación de bosques aledaños a la CDMX afectados por contaminación atmosférica de la convocatoria FORDECYT 2019-5.

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)

Los datos del muestreo corresponden a los datos colectados con kobo y limpiados previamente con el script 1_preprocesamiento_datos_kobo.Rmd.

# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)
parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)

# pivot long parcelas data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy, 
                            cols = healthy:worm, 
                            names_to = "tree_health_simplified",
                            values_to = "n_trees")

Distribución del estado de salud de los árboles por parcela

La siguiente figura muestra el total de árboles muestreados en cada parcela de 10x10 m, y cuántos de estos están bajo alguna categoría de daño:

# Make a nice color pallete and legend order for all plots

my_cols=c("darkgreen", 
              "darkred", 
              "orangered1", 
              "cadetblue", 
              "tan", 
              "beige", 
              "burlywood4", 
              "coral", 
              "aquamarine3", 
              "gray70", 
              "black")

desired_order=c("healthy", 
                "ozone", 
                "ozone_and_other", 
                "others_combined", 
                "drougth", 
                "fungi", 
                "insect", 
                "worm", 
                "acid_rain", 
                "other", 
                "dead")
  
 spanish_labels=c("Sano", 
                  "Ozono", 
                  "Ozono y otros", 
                  "Otros combinados no-ozono", 
                  "Sequía", 
                  "Hongos", 
                  "Insectos", 
                  "Gusano de seda", 
                  "Lluvia acida", 
                  "Otro", 
                  "Muerto")
p <- ggplot(parcelas_long, aes(x=plot, y=n_trees,     fill=tree_health_simplified)) +
  geom_bar(stat="identity") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") 
  

p + theme_bw() +
  ggtitle("Estado de salud de los árboles por parcela") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="Parcelas", y= "número de árboles") +
  theme(text = element_text(size = 25))

Esta es la distribución de las 48 parcelas:

# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

## configure google api

# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).

# if you get the following error when running get_map():

#"Error in aperm.default(map, c(2, 1, 3)) : 
#  invalid first argument, must be an array " 

# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

##  load and register api
api <- readLines("api_key.api")
register_google(key = api)

## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")

# plot sampled plots
p_satmap <-  ggmap(sat_map)
p_satmap + geom_point(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude),
                      color="red") +
          geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                 #    check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5) +
    theme(text = element_text(size = 30))

Según nuestro muestreo, el daño por ozono se distribuye de esta forma:

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = 1.5,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 20))

Haciendo un acercamiento a la zona central:

# zoom plot
sat_map = get_map(location = c(lon = -99.303, lat = 19.2890), zoom = 16, maptype = 'satellite', source = "google")

# plot
p_satmap <-  ggmap(sat_map)

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = .7,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 30)) +
  
  # add plots id
geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                     check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5)

Podemos analizar el daño por ozono en términos de porcentaje de árboles dañados por parcela.

Según la altitud:

# Create new variable with porcentage of ozonoe damage
parcelas_tidy<-parcelas_tidy %>% rowwise() %>% 
                     mutate(., 
                      total=sum(healthy,ozone,ozone_and_other,
                          drougth, acid_rain, other,
                          others_combined, dead, fungi,
                          insect, worm)) %>%
                    mutate(perc.ozone= sum(ozone, ozone_and_other)/total)

#plot
p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_altitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según altitud") + 
  labs(x="Altitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

O según la latitud. Esto es relevante porque latitudes más al sur están más lejos de la CDMX y por ende de la fuente de contaminantes:

p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_latitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según latitud") + 
  labs(x="Latitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

Distribución del estado de salud de árboles individuales

Examinemos el daño por ozono dependiendo de si la planta fue reforestada o no, y de si se encuentra cubierta o expuesta.

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Número de árboles") +
  theme(text = element_text(size = 20))

En términos porcentuales:

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count", position = "fill") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Porcentaje de árboles") +
  theme(text = element_text(size = 20))

Examinemos el daño según el tamaño de los árboles, lo cual es un reflejo de su edad aproximada:

p <- ggplot(muestreo_tidy) +
 scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()
p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de los árboles según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

Explorar con mayor detalle los árboles menores a 15 metros:

p <- filter(muestreo_tidy, tree_heigth<15, tree_nodes>0) %>% 
     ggplot(.) +
     scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()
p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))

Podemos ver lo mismo, pero con el número de nodos:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified))  +
    labs(x="Número de nodos", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

Y si fueron reforestadas o no:

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

En términos porcentuales:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified),
                      position= "fill", binwidth=1)  +
    labs(x="Número de nodos", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 110 rows containing missing values (geom_bar).

También lo podemos explorar por categorías diamétricas:

# changer order of levels of tree_diameter_category
muestreo_tidy$tree_diameter_category<- factor(muestreo_tidy$tree_diameter_category,
levels=c("", "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")) 

levels(muestreo_tidy$tree_diameter_category)<-c(NA, "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count")  +
    labs(x="Categoría diamétrica", y= "Número de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

En términos porcentules:

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count", position = "fill")  +
    labs(x="Categoría diamétrica", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

Ahora evaluaremos el nivel del daño por ozono, es decir, qué porcentaje del árbol se encuentra dañado, dentro del subset de árboles que están dañados.

my_cols2<-c("gold2", "chocolate1", "orangered", "red4", "darkorchid4")
desired_order_percentage<-c("less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")

p_od<- muestreo_tidy %>% filter(!is.na(ozone_damage_percentage)) %>%
            ggplot() +
            scale_fill_manual(values= my_cols2, 
                              breaks = desired_order_percentage,
                              name= "Porcentaje del árbol \n dañado por ozono") +
            theme_bw() + theme(text = element_text(size = 20)) 

p_od +
  geom_bar(aes(x=tree_diameter_category,
               fill=ozone_damage_percentage)) +
  labs(x="Categoría diamétrica", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado en árboles con daño por ozono") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))

Viéndolo por la edad de los árboles (en árboles <15 años):

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage)) +
  labs(x="Número de nodos", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 102 rows containing non-finite values (stat_count).

En términos porcentuales:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position = "fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 102 rows containing non-finite values (stat_count).

Y dividido por reforestadas y naturales:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 102 rows containing non-finite values (stat_count).

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ tree_health_simplified, 
             labeller = as_labeller(c("ozone" = "Sólo daño por ozono",
                                    "ozone_and_other" = "Daño por ozono y otro")))
## Warning: Removed 102 rows containing non-finite values (stat_count).